To do this exercise you will need:
The Carbon file will be in the directory Inputs/Graphene
The Boron and Nitrogen files will be in the directory Inputs/hBN
#Import the libraries we will be needing
import numpy as np
import sisl
import sisl.viz
import plotly.express as px
from pathlib import Path
import os
import itertools
import matplotlib.pyplot as plt
Now we will generate the graphene using sisl
# Create the geometry
graphene = sisl.geom.graphene()
Now we will write this structure in a .fdf file
inputs_dir = Path("Inputs/Graphene")
graphene.write(inputs_dir / "graphen_geom.fdf")
Now we define (and create) the directory where we will be studying the effect on changing the k-grid
# Define the root directory for the convergence study for k grid points
root = Path("k_study")
# Create the directory if it doesn't exist
root.mkdir(exist_ok=True)
Now we create a loop that:
# Set a list with all the k-grid size (in x and y) we want to try
kgrid_points = ['1', '2', '3', '4', '5', '6','7', '8', '9', '12', '15', '18','21', '24']
Es=[] # Where we will store the total energies of each system
tComp=[] # Where we will store the Elapsed time
for kgrid_points in kgrid_points:
# With this part we will create each directory
k_dir = root / kgrid_points
k_dir.mkdir(exist_ok=True)
#print(f"Working on directory with kdrid size of {kgrid_points}") --> for debugging
# With this part we copy the input in each directory
os.system(f"cp {inputs_dir}/* {k_dir}")
# With this part we create the RUN.fdf file with every part
with open(k_dir / 'RUN.fdf', 'w') as f:
f.write("%include graphen_geom.fdf \n") # Includes the geometry file
f.write(f"PAO.BasisSize DZP\n") # Includes the PAO.Basis
#We now include the kgrid
f.write(f"%block kgrid.MonkhorstPack \n")
f.write(f" {kgrid_points} 0 0 0 \n")
f.write(f" 0 {kgrid_points} 0 \n")
f.write(f" 0 0 1 0 \n")
f.write(f"%endblock kgrid.MonkhorstPack")
# Now we run the simulation
siesta_path = "/Users/danibedmarromero/anaconda3/envs/siesta/bin/siesta"
os.system(f"cd {k_dir}; {siesta_path} RUN.fdf > RUN.out")
# Now we store the energy of the system in the Es list
ESystem = sisl.get_sile(k_dir / "RUN.out").read_energy()['total']
Es.append(ESystem)
#print(f"Energy of the system= {ESystem} eV")
# Now we store the Elapsed time in the list tComp
timer =[]
myfile = open(k_dir / 'RUN.out', 'r')
myline = myfile.readline()
while myline:
if 'Elapsed' in myline:
timer.append(myline)
myline = myfile.readline()
myfile.close()
Timer = timer[0]
ElTime = float(Timer[39:44:1])
tComp.append(ElTime)
kgrid_points = ['1', '2', '3', '4', '5', '6','7', '8', '9', '12', '15', '18','21','24']
print('Done!')
Gamma-point calculation with interaction between periodic images Some features might not work optimally: e.g. DM initialization from atomic data Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed Job completed
Done!
Job completed
Now we represent how the energy of each system and the computation time of it change according to the size of the k-grid:
figure = px.line(
x=kgrid_points, y=Es,markers=True
).update_layout(
title="Convergence of energy and its computational cost with k grid points",
xaxis_title="k grid points in the x and y directions",
yaxis_title="Energy [eV]"
)
figure.add_scatter(
x=kgrid_points, y=tComp, yaxis="y2", line_color="red", name=r"Time [s]"
).add_scatter(
x=kgrid_points, y=Es, yaxis="y1", line_color="blue", name=r"Energy [eV]"
).update_layout(
yaxis2 = {"side":'right', 'overlaying': 'y', 'title': "Time [s]"}
).update_layout(
yaxis = {"side":'left', 'title': "Energy [eV]"}
)
figure.show()
We can observe that as the size of the k-grid increases, the total energy evolves towards a constant value (308 eV), which does not change even as the grid size increases. This is a very good indicator that the results have converged because despite increasing our precision, we do not observe considerable changes. In this case, we could say that from a k-grid size of around 7x7 in the x and y directions, we already have enough (even with 4x4 we could be satisfied)!
If we observe the computation time, we notice a very peculiar trend: first, the time decreases until reaching a minimum and then increases again. The second increase in computation time was expected, as an increase in the precision with which we sample the k-space would naturally increase the computational cost. However, the fact that the cost for simulations with low precision is higher than for simulations with an acceptable precision (according to the energy convergence) is surprising. I believe this is due to the system not being able to evolve naturally under the given conditions, leading to an anomalous behavior that is more costly to simulate.
Now we need to simulate the graphene and the hBN to compare them
# Create the hBN by modifying the graphen structure
graphene = sisl.geom.graphene()
hBN = graphene.move([0.03, 0, 0], atoms=1)
#Write the geometry in a .fdf file
input_dir_hBN = Path("Inputs/hBN")
hBN.write(input_dir_hBN / "hBN_geom.fdf")
input_dir_graphene = Path("Inputs/graphene_study")
graphene.write(input_dir_graphene / "graphene_geom.fdf")
Now we need to modify the file, because we currently have this:
And we need to have this:
# Define inputs paths
input_dir_hBN = Path("Inputs/hBN")
# Define the root directory for the hBN study
root_hBN = Path("hBN_study")
# Create the directory if it doesn't exist
root_hBN.mkdir(exist_ok=True)
# Define inputs paths
input_dir_graph = Path("Inputs/graphene_study")
# Define the root directory for the hBN study
root_graph = Path("Graphene_study")
# Create the directory if it doesn't exist
root_graph.mkdir(exist_ok=True)
#Copy the hBN input to where we will run the simulation
os.system(f"cp {input_dir_hBN}/* {root_hBN}")
#Creamos el fichero Run.fdf
with open(root_hBN / 'RUN.fdf', 'w') as f:
# We include the fdf file that contains the geometry, with the %include statement
# Note the \n, which means "new line" in text files.
f.write("%include hBN_geom.fdf \n")
# We now include the basis size input.
f.write(f"PAO.BasisSize DZP\n")
#We now include the kgrid
f.write(f"%block kgrid.MonkhorstPack \n")
f.write(f" 7 0 0 0 \n")
f.write(f" 0 7 0 0 \n")
f.write(f" 0 0 1 0 \n")
f.write(f"%endblock kgrid.MonkhorstPack\n")
f.write("SaveRho true\n")
f.write("TS.HS.Save true\n")
#Copy the hBN input to where we will run the simulation
os.system(f"cp {input_dir_graph}/* {root_graph}")
#Creamos el fichero Run.fdf
with open(root_graph / 'RUN.fdf', 'w') as f:
# We include the fdf file that contains the geometry, with the %include statement
# Note the \n, which means "new line" in text files.
f.write("%include graphene_geom.fdf \n")
# We now include the basis size input.
f.write(f"PAO.BasisSize DZP\n")
#We now include the kgrid
f.write(f"%block kgrid.MonkhorstPack \n")
f.write(f" 7 0 0 0 \n")
f.write(f" 0 7 0 0 \n")
f.write(f" 0 0 1 0 \n")
f.write(f"%endblock kgrid.MonkhorstPack\n")
f.write("SaveRho true\n")
f.write("TS.HS.Save true\n")
#Run the simulations
siesta_path = "/Users/danibedmarromero/anaconda3/envs/siesta/bin/siesta"
os.system(f"cd {root_hBN}; {siesta_path} RUN.fdf > RUN.out")
Job completed
0
#Run the simulations
siesta_path = "/Users/danibedmarromero/anaconda3/envs/siesta/bin/siesta"
os.system(f"cd {root_graph}; {siesta_path} RUN.fdf > RUN.out")
Job completed
0
Proceed with the analysis
RUN_hBN = Path('hBN_study/RUN.FDF')
H_hBN = sisl.get_sile(RUN_hBN).read_hamiltonian()
RUN_graph = Path('Graphene_study/RUN.FDF')
H_graph = sisl.get_sile(RUN_graph).read_hamiltonian()
# We need to define a path of k points
band_struct_hBN = sisl.BandStructure(H_hBN, points=[[0, 0, 0], [2/3, 1/3, 0], [1/2, 0, 0]],
divisions=100, names=[r"\Gamma", "M", "K"]
)
# Then we can plot the bands
band_struct_hBN.plot()
# We need to define a path of k points
band_struct_graph = sisl.BandStructure(H_graph, points=[[0, 0, 0], [2/3, 1/3, 0], [1/2, 0, 0]],
divisions=100, names=[r"\Gamma", "M", "K"]
)
# Then we can plot the bands
band_struct_graph.plot()
With the band structure plots, we can appreciate that the main difference (or at least with more direct practical effects on the compound) is that in hBN we find a GAP (of about 4.5 eV), while in graphene we find its famous Dirac cone, which grants it its very special behavior.
# Get the fatbands plot
fatbands_hBN = band_struct_hBN.plot.fatbands()
# Split the contributions by the n and l quantum numbers
fatbands_hBN.split_groups(on="n+l")
# Get the fatbands plot
fatbands_graph = band_struct_graph.plot.fatbands()
# Split the contributions by the n and l quantum numbers
fatbands_graph.split_groups(on="n+l")
With the fatbands, we can observe how the orbitals n=2, l=1 and n=3, l=2, which give rise to the Dirac cone (especially the orbital n=2, l=1), look like. These same orbitals are the ones that contribute to the bands forming the GAP in hBN.
# Get the PDOS plot
pdos_plot_hBN = H_hBN.plot.pdos(
kgrid=[90,90,1], nE=1000, Erange=[-10, 10],
distribution={"method": "gaussian", "smearing": 0.1}
)
# Split the contributions by the n and l quantum numbers
pdos_plot_hBN.split_DOS(on="n+l", name="Atom $atoms")
# Get the PDOS plot
pdos_plot_graph = H_graph.plot.pdos(
kgrid=[90,90,1], nE=1000, Erange=[-10, 10],
distribution={"method": "gaussian", "smearing": 0.1}
)
# Split the contributions by the n and l quantum numbers
pdos_plot_graph.split_DOS(on="n+l", name="Atom $atoms")
With the pDOS plots, it's very easy to see how in graphene we have orbitals (the p) that contribute to the density of states, whereas in hBN we have a GAP.
rho_hBN = sisl.get_sile(RUN_hBN).read_grid("RHO")
rho_hBN.plot(axes="xy", plot_geom=True)
rho_graph = sisl.get_sile(RUN_graph).read_grid("RHO")
rho_graph.plot(axes="xy", plot_geom=True)
We can appreciate how the electron density is distributed differently, as even though the centers of the hexagons are where the most negative charge density is located, in graphene we have much more charge density on the bonds (thanks to the pi bonds), while in hBN the charge density is completely localized on the Nitrogen atoms. This is because the bonds have a strong donor-acceptor character, as Nitrogen is much more electronegative than Boron.